Minimization of an energy error functional to solve 
a Cauchy problem arising in plasma physics: the 
reconstruction of the magnetic flux in the vacuum 
surrounding the plasma in a Tokamak 
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Abstract. A numerical method for the computation of the magnetic flux in the 
vacuum surrounding the plasma in a Tokamak is investigated. It is based on the 
formulation of a Cauchy problem which is solved through the minimization of an 
energy error functional. Several numerical experiments are conducted which show the 
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1. Introduction 

In order to be able to control the plasma during a fusion experiment in a Tokamak it is 
mandatory to know its position in the vacuum vessel. This latter is deduced from the 
knowledge of the poloidal flux which itself relies on measurements of the magnetic field. 
In this paper we investigate a numerical method for the computation of the poloidal 
flux in the vacuum. Let us first briefly recall the equations modelizing the equilibrium 
of a plasma in a Tokamak [1] . 

Assuming an axisymmetric configuration one considers a 2D poloidal cross section 
of the vacuum vessel Qy in the (r,z) system of coordinates (Fig. [1]). In this 

setting the poloidal flux tp(r, z) is related to the magnetic field through the relation 

1 dtp dtp 
(B r ,B z ) = -(-t-itt) an d> as there is no toroidal current density in the vacuum 

r dz or 
outside the plasma, satisfies the following equation 

Ltp = in Q x (1) 

d 1 d. d 1 d. 
where L denotes the elliptic operator L. = — \tt(~tt) + "^-( - t^)1 and Qx = ^V — ^p 

or r or oz r oz 
denotes the vacuum region surrounding the domain of the plasma Qp of boundary Tp 

(see Fig. [2]). Inside the plasma Eq. ([1]) is not valid anymore and the poloidal flux 

satisfies the Grad-Shafranov equation J2J [3] which describes the equilibrium of a plasma 

confined by a magnetic field 

Ltp = fi j(r,tp) in Q P (2) 

where /j,q is the permeability of the vacuum and j(r, tp) is the unknown toroidal current 
density function inside the plasma. Since the plasma boundary Tp is unknown the 
equilibrium of a plasma in a Tokamak is a free boundary problem described by a 
particular non-linearity of the model. The boundary is an iso-flux line determined 
either as being a magnetic separatrix (hyperbolic line with an X-point as on the left 
hand side of Fig. [2]) or by the contact with a limiter (Fig. [2] right hand side). In other 
words the plasma boundary is determined from the equation tp(r, z) = tpp, tpp being 
the value of the flux at the X-point or the value of the flux for the outermost flux line 
inside a limiter. 

In order to compute an approximation of tp in the vacuum and to find the plasma 
boundary without knowing the current density j in the plasma and thus without using 
the Grad-Shafranov equation ([2]) the strategy which is routinely used in operational 
codes mainly consists in choosing an apriori expansion method for tp such as for example 
truncated Taylor and Fourier expansions for the code Apolo on the Tokamak ToreSupra 
[I] or piecewise polynomial expansions for the code Xloc on the Tokamak JET [5l E]. 
The flux tp can also be expanded in toroidal harmonics involving Legendre functions or 
expressed using Green functions in the filament method ([HE], [9] and the references 
therein). In all cases the coefficients of the expansion are then computed through a fit 
to the measurements of the magnetic field. Indeed several coils and flux loops surround 
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Figure 1. Cross section of the vacuum vessel. The domain ily, its boundary Ty. 
Coils providing measurements of the components of the magnetic field tangent and 
normal to Ty are represented surrounding the vacuum vessel. 



limiteF 




Figure 2. The plasma domain fip and the vacuum region VLx- The plasma boundary 
is determined by an X-point configuration (left) or a limiter configuration (right). The 
fictitious contour T/ is represented inside the plasma. 



the boundary Ty of the vacuum vessel and measure the magnetic field and the flux (see 

Fig. m>. 

In this paper we investigate a numerical method based on the resolution of a Cauchy 
problem introduced in ([ID], Chapter 5) which we recall here below. The proposed 
approach uses the fact that after a preprocessing of these measurements (interpolation 
and possibly integration on a contour) one can have access to a complete set of Cauchy 

data, / = -0 on T v and g = ——— on T v . 

r on 



A Cauchy problem in plasma physics 

The poloidal flux satisfies 

Lip = in fix 

on rv 



(3) 



g on r 



v 



4> = f 

I dip _ 

r dn 

ip = ipp on Tp 

In this formulation the domain fix — ^x(^) is unknown since the free plasma 
boundary Tp as well as the flux ipp on the boundary are unknown. Moreover the problem 
is ill-posed in the sense of Hadamard [11] since there are two Cauchy conditions on the 
boundary T V - 

In order to compute the flux in the vacuum and to find the plasma boundary we 
are going to define a new problem as in [10] which is an approximation of the original 
one. Let us define a fictitious boundary Ti fixed inside the plasma (see Fig. [2]). We 
are going to seek an approximation of the poloidal flux tp satisfying Lip = in the 
domain contained between the fixed boundaries IV and Tj. The problem becomes one 
formulated on a fixed domain fl: 
( Lip = in Q 



ip = f on I\ 



(4) 



I dip 
r dn 



on r 



v 



However no boundary condition is known on Tj. One way to deal with this second 
issue and to solve such a problem is to formulate it as an optimal control one. Only the 
Dirichlet condition on iy is retained to solve the boundary value problem and a least 
square error functional measuring the distance between measured and computed normal 
derivative and depending on the unknown boundary condition on T/ is minimized. Due 
to the illposedness of the considered Cauchy problem a regularization term is needed 
to avoid erratic behaviour on the boundary where the data is missing. A drawback 
of this method developed in [10] is that Dirichlet and Neumann boundary conditions 
on ry are not used in a symmetric way. One is used as a boundary condition for the 
partial differential equation, Lip = 0, whereas the other is used in the functional to be 
minimized. 

It should be noticed here that by the Cauchy-Kowalewska theorem [TT] the function 
ip can be extended in the sense of Lip = in a neighborhood of Tp inside the plasma. 
Hence the problem formulated on a fixed domain with a fictitious boundary Ti not 
"too far" from Tp is an approximation of the free boundary problem. As mentioned 

in [TO] if Ti were identical with Tp then by the virtual shell principle [12] the quantity 

ldip t . 1 . 

w = -Ti - rv would represent the surface current density (up to a factor — ) on Tp for 
r on iiq 

which the magnetic field created outside the plasma by the current sheet is identical to 
the field created by the real current density spread throughout the plasma. 
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Freezing the domain to Q by introducing the fictitious boundary Tj enables to 
remove the nonlinearity of the problem. The plasma boundary Tp can still be computed 
as an iso-flux line and thus is an output of our computations. We are going to compute 
a function ip such that the Dirichlet boundary condition u = -0 on Tj is such that the 
Cauchy conditions on IV are satisfied as nearly as possible in the sense of the error 
functional defined in the next Section. 

The originality of the approach proposed in this paper relies on the use of an error 
functional having a physical meaning: an energy error functional or constitutive law 
error functional. Up to our knowledge this misfit functional has been introduced in [13] 
in the context of aposteriori estimator in the finite element method. In this context, the 
minimization of the constitutive law error functional allows to detect the reliability of the 
mesh without knowing the exact solution. Within the inverse problem communauty this 
functional has been introduced in [HI [151 03] in the context of parameter identification. 
It has been widely exploited in the same context in |17j . It has also been used for 
Robin type boundary condition recovering [18] and in the context of geometrical flaws 
identification (see [19] and references therein). For lacking boundary data recovering 
(i.e. Cauchy problem resolution) in the context of Laplace operator, the energy error 
functional has been introduced in [20| |2"T] . A study of similar techniques can be found 
in [22l [23] and the analysis found in these papers uses elements taken from the domain 
decomposition framework [24|. 

The paper is organized as follows. In Section 2 we give the formulation of the 
problem we are interested in and provide an analysis of its well posedness. Section 3 
describes the numerical method used. Several numerical experiments are conducted to 
validate it. The final experiment shows the reconstruction of the poloidal flux and the 
localization of the plasma boundary for an ITER configuration. 

2. Formulation and analysis of the method 

2.1. Problem formulation 

As described in the Introduction the starting point is the free boundary problem ([3]). 
We first proceed as in [10] and in a first step consider the fictitious contour Tj fixed 
in the plasma and the fixed domain Q contained between Ty and Tj. Problem ([3]) is 
approximated by the Cauchy problem (jlj). 

In a second step the problem is separated into two different ones. In the first one 
we retain the Dirichlet boundary condition on Ty only, assume v is given on Tj and 
seek the solution ip D of the well-posed boundary value problem: 

LipD = in Q 

< *Pd = f on Ty (5) 

k i\) D = v onT/ 
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The solution tp D can be decomposed in a part linearly depending on v and a part 
depending on / only. We have the following decomposition: 

fa = M", f) = M", o) + <Mo, /) := Mv) + Mf) (6) 

where 4>d( v ) an d 4>D(f) satisfy: 
Lifjrj{v) = in fl 



i> D {v) 



onf. 



f i4 D (f) = infi 

"ipD(f) = f on T y 
I Mf) = on r / 



(7) 



k iPd(v) = v on Tj 

In the second problem we retain the Neumann boundary condition only and look 
for -07V satisfying the well-posed boundary value problem: 

LipN = in Q 

ldip N 



g on T 



v 



(8) 



r dn 

■ipN = v on T/ 

in which ijj^ can be decomposed in a part linearly depending on v and a part 
depending on g only We have the following decomposition: 

ipN = 4>n{v, g) = iPn{v, 0) + ip N (0, g) := ip N (v) + ^ N {g) (9) 



where 



Lipj^iv) = in Q 

±dip N (v) 

= on Ti 



Lifj N (g) =0 in Q 

1 9 ^ r f io) 

g> on 1 v ^ ' 



r dn 
^ tp N = on T/ 



r <9n 

i/jn(v) = v on T/ 

In order to solve problem (JIJ), / G iJ 1//2 (IV) and (7 G H'~" 1 ' 2 (rV) being given, we 
would like to find u EU = H 1//2 (T/) such that -0 = iPd(u, f) = i^n{u, g). To achieve this 
we are in fact going to seek u such that J(u) = inf J(y) where J is the error functional 

defined by 

J(n)= 1 - f-\\^ D (uJ)-Vij N (u } g)\\ 2 dx (11) 

1 Jn r 

measuring a misfit between the Dirichlet solution and the Neumann solution. 



2.2. Analysis of the method 

In order to minimize J one can compute its derivative and express the first order 
optimality condition. When doing so the two symmetric bilinear forms sd and sn 
as well as the linear form I defined below appear naturally and in a first step it is 
convenient to give a new expression of functional (ITT1) using these forms. 
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Let u,v e // 1/2 (r 7 ) and define 



s D (u,v)= / -Vi/j D (u)V^ D (v)dx (12) 

Jnr 

Applying Green's formula and noticing that iPd{v) = v on Tj and iI>d(v) = on Ty we 
obtain 

s D (u,v)= / -d n if; D (u)i/j D (v)do-- V(~Vij} D (u))i/) D (v)dx = I -d n ip D (u)vda(13) 
Jan r Jn r J Tl r 

where the integrals on the boundary are to be understood as duality pairings. In Eq. f|T3|) 
one can replace iPd(v) by any extension TZ(v) in Hq(Q, T v ) = {^ 6 H 1 (Q),ip^ v = 0} of 
v G H^iYj). 
Hence sp can be represented by 

s D (u,v) = / -Vif> D (u)VK(v)dx (14) 

Jnr 

Equivalently sn is defined by 

s N (u,v) = / -VifjNMVi/jNMdx (15) 

Jnr 

Since ipisriv) = v on Ti and -dnijj^iu) = on Ty we have that 

r 

Sn(u, v) — I -d n lfj N (u)lp N (v)d<J— V(-VV ; Ar(M))V ; Ar(f )dx = / -(9„'i/ ; Af(w)^^Cr(16) 

J an r Jn r Vrv r 

and sn can also be represented by 

s N (u,v)= -Vi/j N {u)VK(v)dx (17) 



where TZ(v) is any extension in H x (f2) of f G Z/ 1//2 (r/). 
Let us now introduce 

F( Wj v) = I f -(V^ D (u,f)-V^ N (u,g))(VMvJ)-Vij N (v,g))dx(l8) 



2 Jnr 
such that J(v) = F(v,v) and the linear form I defined by 

Kv) = ~ ( k^Mf) - V*p N (g))VMv)dx (19) 

Jnr 

which can also be computed as 

Kv) = - [ -Wv(f) ~ Vi> N (g))VK(v)dx (20) 

Jnr 

Using the decompositions of Eqs. ((6]) and (|9]) as well as the fact that for all 
v G H^iTj) 

-VMf)Vi/j N (v)dx = [ Mf)-d n i/} N (v))da- [ $ D (f)V(-ViJj N (v))dx = 0(21) 
n r Jon r Jn r 



and 



~V$ N {g)V^ N {v)dx = (22) 

nr 
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and finally that for all u and v G H l l 2 (Yi) 



-Vip N {u)Vilj N (v) = I -Vip N {u)Vij D {v) = [ -Vij D (u)V<iP N (v) (23) 
n r Jn r J n r 

it is shown with that 

F{u, v) = -{s D {u, v) - s N (u, v) - l(u) - /(f)) + c (24) 

where the constant c is given by 

c=\ [-\\VMf)-^ N (9)\\ 2 dx (25) 

1 Jn r 

Hence functional J can be rewritten as 

J{ v ) = 7)( s d{v, v) - s N (v, v)) - l(v) + c (26) 

Following the analysis provided in [22] it can be proved that in the favorable case 
of compatible Cauchy data (/, g) the Cauchy problem admits a solution. There exists 
a unique u 6 U such that ip D (u, f) = ipN(u,g). The minimum of J is also uniquely 
reached at this point, J(u) = 0. This solution is given by the first order optimality 
condition which reads 

(J'(u),v) = s D {u,v)-s N {u,v) -l{v) = \fveU (27) 

Equation (|27|) has an interpretation in terms of the normal derivative of ipr> and iftpj on 
the boundary From Eqs. (TT3|) and (!T6|) and from 

l(v) = - f -(V#D(f)-VJ> N (g))VMv)dx = - / -(d n Mf)-d n $ N (g))vda(2l 
Jn r J Tl r 

we deduce that the optimality condition can be rewritten as 

/ [{-d n M% f) ~ ~d n Mu, g))\vda = \/veU (29) 

J Tl r r 

which can be understood as the equality of the normal derivatives on Tj. 

This can also be linked with domain decomposition methods and tools: the first 
optimality condition when minimizing J amounts to solve an interfacial equation 

(S D - S N )(v) = x, 

where Sn and S^ are the Dirichlet-to-Neumann operators associated to the bilinear 
forms and defined by: 

S D : H^iTj) — ► H-^iTj) 

IdMv) (30) 



v 



r dn 



s N : F 1 /2 (r/) _^ ff-i/a(r 7 ) 

lfy N (v) (31) 



v 



r dn 



\di) D ld^ N 

and x = ^ ! TT~ on r ^- 

r on r on 
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When dealing with domain decomposition tools to solve such a PDE, one has to 
solve the interfacial equation on Tj (see an illustration of domain decomposition on Fig. 

ED 

S(v) = x, 

where S := Sd + Sn (the sign + is related to the change of unit normal orientation). 
On the other hand in the missing boundary data recovering context one has to solve 



S(v) 



X- 



with S :— Sd — Sn- In this case Qi and f2 2 are equal to the whole domain Q and we 
resort here to a fictitious domain decomposition process (see Fig. [3]). 



A i 



II 



Q 1 



% 



a, 




Figure 3. Illustration of a domain decomposition method and of a fictitious domain 
decomposition method for data completion 



Since Sd and Sn have the same eigenvectors and have asymptotically the same 
eigenvalues, the interfacial operator S = Sd — Sn is almost singular [22]. This point 
together with the fact that the set of incompatible Cauchy data is known to be dense in 
the set of compatible data (and thus numerical Cauchy data can hardly by compatible) 
make this inverse problem severely ill-posed. 

Some regularization process has to be used. One way to regularize the problem 
is to directly deal with the resolution of the underlying quasi- singular linear system 
using for example a relaxed gradient method [201 [2T] . In this paper we have chosen a 
regularization method of the Tikhonov type. It consists in shifting the spectrum of S 
by adding a term 

(Sd — Sn) + sSd- 

where e is a small regularization parameter. This regularization method is quite natural 
since the ill-posedness of the inverse problem and the lack of stability in the identification 
of u by the minimization of J is strongly linked to the fact that J is not coercive (see 
[22] and below). We are thus going to minimize the regularized cost function: 



with 



J £ (v) = J(v) + eR D (v) 

R D (v) = l[-\\VMv)\\ 2 dx 
This brings us to the framework described in [23]. We want to solve the following 
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Problem P £ : find u £ EU such that J £ (u £ ) = inf JJv) 

and the following result holds. 

Proposition 1 (i) Problem P £ admits a unique solution u £ G W characterized by the 
first order optimality condition 

(J' £ (u £ ),v) = es D (u £ ,v) + s D (u £ ,v)-s N (u £ ,v)-l(v) = \/v G W(32) 

(ii) For a fixed e the solution is stable with respect to the data f and g. 
Iff 1 , P e H l l 2 (T v ) and g 1 , g 2 G H- l / 2 {V v ) it holds that 

C 

\\ u l ~ ulWrn/'^Tj) < jdlf ~ f 2 \\m/Hr v ) + Wg 1 - 9 2 \\h-^(t v )) (33) 

(Hi) If there exists u eW such that iPd{ u , f) — ^jv(w, g) then u £ — >■ u inU when e — > 0. 



proof: 
1. We need to prove the continuity and the coercivity of J £ . 
Continuity. 

The maps v h- y i/jd(v) and v i-> iPn{v) are continuous and linear from H 1 ^ 2 (Ti) to H 1 ^). 
Moreover since "0d(/) and i^N{g) are in H l (Q) and Tm > r > r m > in f2 it is shown 
with Cauchy Schwartz that the bilinear forms s^ and sn, the linear form I and thus J £ 
are continuous on H l l 2 (Ti). 

Coercivity. 
The bilinear form so is coercive on ^^(Tj). Indeed 4>d{v) G i?o(f2, T^) an d 
the Poincare inequality holds. There exists en > such that c^HV^i^i;)^/™ > 
\\Mv)\\h {n y Thus 



s D (v,v)= / -||V« V )|| 2 >— / ||V^(t;)ir> d l^frOllfngi) 

Jn ^ r M Jn rM\i- ~r en) 

The application ipo{v) G H l {Vt) — > 4>D(v)\rj = v G i/ 1 / 2 ^/) is continuous and 

lkll^i/2 (ri) < c'nllV'D^lllim). Finally we obtain s D (w,w) > - J HIgi/^r,)- 

On the contrary, since for iJjn(v) G i? 1 (fi) the seminorm does not bound the L 2 norm, 
the bilinear form s^ is not coercive and because of the minus sign in s = sp — sn we 
need to prove that s(v, v) > to obtain the coercivity of the bilinear part of functional 
J £ . We use the same type of argument as in [22] : 

Let v G H 1 ^) be such that u|r 7 = v. Define the symmetric bilinear form 

a(i/j,(/)) = / -VipVifidx. The decomposition ipn(v) = tp% + v holds. ip% is the unique 



solution of the variational formulation: 

find ip° D G H£(n) such that a(ip%, </>) = a(v, 0) V0 G H^(Q). 

Since the bilinear form a is symmetric ^ is also the solution to a minimization 
problem, 

ip D = argmin V(0) 
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where V{4>) = -a((j), 4>) + a(v, (p). 

In the same way it holds that iJ)n(v) — 4>n + ^ with 

ij) N = argmin V((f>) 
0eHj(fi,r/) 

Since i/j° d G H%(Q) C £#(fi,Ij) then V(ip%) < V(ip° D ). Thus 

a(V&, ip° N ) + 2a(v, ip° N ) + a(v, v) < a(i/j° D , i/j° d ) + 2a(v, tjP D ) + a(v, v) 

i.c 

a{il> N {v),il) N (v)) < a{il) D {v),ij) D {v)) 

and 

sn(v,v) < s D (v,v) 

Eventually it holds that 

-s(v,v) + -s D (v,v) > Ce\\v\\ 2 Hl/2{Vi) 

Using the continuity and the coercivity of J £ it results from [25] that problem P £ admits 
a unique solution u £ G U. 

The solution u £ is characterized by the first order optimality condition which is 
written as the following well-posed variational problem 

(J f e (u e ),v) = es D (u £ ,v) + s D (u £ ,v) - s N (u £ ,v) - l(v) = \fv G U (34) 

which as in Eq. ( )29|) can be understood as an equality on Yj. 



2. The stability result is deduced from the optimality condition (pK 

Let u\ (resp. u 2 ) be the solution associated to (Z 1 ,*? 1 ) (resp. (f 2 ,g 2 )). The 
optimality condition reads 

esn(ul,v) + sn(ul,v) — sx(ul,v) — h(v) = Q Vv <E U (35) 

es D (u 2 , v) + s D (u 2 , v) — s N (u 2 ,v) — l 2 (v) = VwGW (36) 

Substracting Eq. (l36l) from ( J35l) . choosing v = u] — u 2 and using the coercivity leads to 

Cell^ - u 2 £ \\ 2 H i/2 (ri) < \(h - h){u] - u 2 )\ 

The map / i-> ipo(f) is linear and continuous from if 1 ' 2 (IV) to if 1 (fi), and so is the 
map g h-> ip^ig) from i7 _1 / 2 (IV) to H 1 ^). Using these facts and Cauchy Schwarz it 
follows that 

l (/l _ y (u i _ tt 2)| < — / Kvfec/ 1 - / 2 ) - v^^ 1 - /)) 

S/ipoiu). — u 2 )\alx 

c 

< — (ll/ 1 - / HjyV^rv) + \\g - 9 2 \\H-w*{r v )) 

I'm 

111 211 

\\ u e ~ U e\\H^{T I ) 
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Eventually 

C 1 
\\ u l -«ellirv<«crj) < ITT^tGI/ 1 ~ f\\H^{v v ) + ll^ 1 - 9 2 \\h-^(t v )) 

3. We will not prove rigourously this last point (the proof of Proposition 3.2 in [23] 
can be adpated). A sketch of the proof is as follows. Let us suppose that there exists 
u &U such that ijJd{u, f) = iJjn{ u , q)- A key point is to show that sd(u £ , u £ ) —¥ Sd(u, u) 
when e —¥ 0. Then in a second step using the optimality conditions for u and u £ it is 
shown that 

s D (u £ - u,u e — u) < s D (u, u) - s D (u £ , u £ ) 

which gives the result thanks to the coercivity of sd in H l l 2 (Ti). 

D 

3. Numerical method and experiments 

3.1. Finite element discretization 

The resolution of the boundary value problems (j^J) and ( TTUj) is based on a classical P 1 
finite element method [26J. 

Let us consider the family of triangulation r^ of Q, and Vh the finite dimensional 
subspace of iJ 1 (fi) defined by 

V h = {tp h G H\ti)^ h \ T G P\T), VT G r h }. 

Let us also introduce the finite element space on Fj 

D h = {v h = ip h \r n i>h e Vh}- 

Consider (0i)j=i,...jv a basis of V^ and assume that the first Ny t mesh nodes (and basis 
functions) correspond to the ones situated on Tj. A function iph G Vh is decomposed as 

i>h = Z)t=l a *^ alld its traCe OI1 r ^ aS v h = i>h\T I = Z)i=l a J0i|lV 

Given boundary conditions t^ on I" 1 / and //,, g^ on Ty one can compute the 
approximations ip D ,h(v h ), ip N ,h(v h ), 4> D ,h(fh) and ^ N ,h(9h) with the finite element 
method. 

In order to minimize the discrete regularized error functional, J £y h{uh) we have to 
solve the discrete optimality condition which reads 

es D>h (u h , v h ) + s Dth (u h , v h ) - s Njh (u h , v h ) - l(v h ) = \/v h G D h (37) 

which is equivalent to look for the vector u solution to the linear system 

Su = 1 (38) 

where the N ri x N ri matrix S representing the bilinear form Sh = ssr>,h + s d,h — SN,h is 
defined by 

Sij = s h ((pi,(f)j) (39) 
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and 1 is the vector {lh{4 > i))i=i,...N v ■ 

In order to lighten the computations the matrices are evaluated by 



sdm 



>j) 



and 



SN,h{(pi,<Pj 



~Vip D>h {<j)i)VK{(j)j)dx 



i)VR{(j)j)dx 



(40) 



where TZ{<f>j) is the trivial extension which coincides with 
elsewhere. 

In the same way the second member 1 is evaluated by 



//. 



-WdAU) - V4>N,h(gh))V1l(<f>i)dx 



(41) 
on Tj and vanishes 

(42) 



It should be noticed here that matrix S depends on the geometry of the problem 
only and not on the input Cauchy data. Therefore it can be computed once for all 
(as well as its LU decomposition for exemple if this is the method used to invert the 
system) and be used for the resolution of successive problems with varying input data 
as it is the case during a plasma shot in a Tokamak. Only the second member 1 has to 
be recomputed. This enables very fast computation times. 

All the numerical results presented in the remaining part of this paper were obtained 



using the software FreeFem++ (http://www.freefem.org/ff++/). We are concerned 
with the geometry of ITER and the mesh used for the computations is shown on Fig. HJ 
It is composed of 1804 triangles and 977 nodes 150 of which are boundary nodes divided 
into 120 nodes on Ty and 30= iV r/ on Yj. The shape of Yj is chosen empirically. 




Figure 4. The mesh used for the ITER configuration in FreeFem++ 



3.2. Twin experiments 

Numerical experiments with simulated input Cauchy data are conducted in order to 
validate the algorithm. Assume we are provided with a Neumann boundary condition 
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function g on T v . We generate the associated Dirichlet function / on T v assuming a 
reference Dirichlet function u re f is known on Tj. We thus solve the following boundary 
value problem: 

Lip Njref (u re f,g) = in fi 

-d n 4> N ,ref(u ref ,g) = g OIlTy (43) 

^N,ref(u re f,g) = u ref on T 7 

and set / = ifj Ntre f(u ref ,g)\r v . 

We have considered two test cases. In the first one (TCI) 

u re f{r, z) = 50 sin(r) 2 + 50 on T 7 (44) 

and in the second one (TC2) u re f is simply a constant 

u re f(r,z) = 40 on T 7 (45) 

The numerical experiments consist in minimizing the regularized error functional 
J £ defined thanks to / and g. The obtained optimal solution u opt and the associated 
ijj opt are then compared to u re f and ip re f which should ideally be recovered. Three cases 
are considered: the noise free case, a 1% noise on / and g and a 5% noise. 

When the noise on / and g is small and the recovery of u is excellent there is very 
little difference between the Dirichlet solution ipD{u pt, f) and the Neumann solution 
^N{uopt, g)- However this not the case any longer when the level of noise increases. The 
Dirichlet solution is much more sensitive to noise on / than the Neumann solution is 
sensitive to noise on g. Therefore the optimal solution is chosen to be ip opt = ipN(u opt , g). 

The results are shown on Figs. [5] and [6] where the reference and recovered solutions 
are shown for the three levels of noise considered. The results are excellent for the noise 
free case in which the Dirichlet boundary condition u is almost perfecty recovered (Fig. 
E|). The differences between u opt and u re f increase with the level of noise (Fig. [7] and 
Tab. CD). As it is often the case in this type of inverse problems the most important 
errors on ip opt are localized close to the boundary Tj and vanishes as we move away from 
it (Fig. ED. 

Tables[2]and[3]sumarize the evolution of the values of J, Rd and J e for the different 
noise level. First guess values (u = 0) are also provided for comparison. Please note 
that the regularization parameter was chosen differently from one experiment to another 
depending on the noise level. This was tuned by hand. In the next section we propose 
to use the L-curve method [27] to choose the value of e. 



3.3. An ITER equilibrium 

In this last numerical experiment we consider a 'real' ITER case. Measurements of 
the magnetic field are provided by the plasma equilibrium code CEDRES++ [28J. 
These mesurements are interpolated (and integrated) to provide / and g on T V - The 
regularized error functional is then minimized to compute the optimal u op t- The choice 
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reference solution 



oplimal solution 0% 





optimal solution 1% 



optimal solution 5% 





Figure 5. First test case (TCI), u re f given by Eq. (|44]l . Top left: reference solution 
^N,ref{u re fTg). Top right: recovered solution with no noise on the data. Bottom left: 
recovered solution with a 1% noise on the data. Bottom right: recovered solution with 
a 5% noise. 



noise level 


error TCI error TC2 


0% 
1% 
5% 


0.0131 
0.0659 
0.1526 


0.0055 
0.0170 
0.0405 



Table 1. Maximum relative error 



^opt 



l ref\ 



for TC 1 and 2 



l ref\ 
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reference solution 



oplimal solution 0% 



optimal solution 1% 




optimal solution 5' 




Figure 6. Second test case (TC2), u re f given by Eq. (J45I) . Top left: reference solution 
ipN,ref(u re f,g). Top right: recovered solution with no noise on the data. Bottom left: 
recovered solution with a 1% noise on the data. Bottom right: recovered solution with 
a 5% noise. 
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Figure 7. u re / and the recovered u op t for the 3 levels of noise on the data. Left: TCI. 
Right TC2. 



A Cauchy problem in plasma physics 



17 




Figure 8. Relative error \ip opt — VVptl/lVVe/l f° r TCI with 5% noise. 





J 


Rd 


Je 


e | 


« = 0no noise 


46.8643 





46.8643 




u opt no noise 


0.0021 


46.8722 


0.0026 


10~ 5 


Uopt 1% noise 


1.8443 


46.5553 


1.8676 


5 x 10" 4 


Uopt 5% noise 


9.2180 


46.5575 


9.2646 


10~ 3 



Table 2. TCI results. Values of the error functional, the regularization term, the 
total cost function and the chosen regularization parameter for the initial guess (row 
1), the optimal solutions for different noise levels (row 2, 3 and 4). 





J 


Rd 


Je 


e | 


« = 0no noise 


30.7231 





30.7231 




u opt no noise 


0.0003 


30.7242 


0.0006 


10~ 5 


Uop t 1% noise 


0.7300 


30.7159 


0.7607 


10~ 3 


u op t 5% noise 


3.6516 


30.6822 


3.8050 


5 x 10- 3 



Table 3. TC2 results. Values of the error functional, the regularization term, the 
total cost function and the chosen regularization parameter for the initial guess (row 
1), the optimal solutions for different noise levels (row 2, 3 and 4). 



of the regularization parameter e is made thanks to the computation of the L-curve 
shown on Fig. [9l It is a plot of (J(u op t)(e),RD(u pt)(£)) as e varies. The corner of the 
L-shaped curve provides a value of e = 5.10~ 4 . 

The computed u opt is shown on Fig. [10] and numerical values are given in Tab. HJ 
The recovered poloidal flux ip is shown on Fig. [TTJ The boundary of the plasma is found 
to be the isoflux if) = 16.3 which shows an X-point configuration. 
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0.0005 '^-0.001 



Figure 9. L-curve computed for the ITER case. The corner is located at e = 5 x 10 




5 



Figure 10. Optimal u op t for the ITER case. 



optimal solution 




Figure 11. Optimal solution for the ITER case. The plasma is in an X-point 
configuration 
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u = 

"'opt 



J Rd J £ 



31.1026 31.1026 

0.8053 39.9169 0.8253 5 x 10" 4 



Table 4. ITER case results. Values of the error functional, the regularization term, 
the total cost function and the chosen regularization parameter for the initial guess 
(row 1) and the optimal solution (row 2) 



4. Conclusion 

We have presented a numerical method for the computation of the poloidal flux in the 
vacuum region surrounding the plasma in a Tokamak. This computation enables in a 
second step the reconstruction of the plasma boundary. The algorithm is based on the 
optimization of a regularized error functional. 

Numerical experiments have been conducted. They show that the method is precise 
and robust to noise on the Cauchy input data. It is fast since the optimization reduces to 
the resolution of a linear system of very reasonable dimension. Successive equilibrium 
reconstructions can be conducted very rapidly since the matrix of this linear system 
can be completely precomputed and only the second member has to be updated. The 
L-curve method proved to be efficient to specify the regularization parameter. 
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